
#include "fe_odens.h"

#include "glodef.h"
#include "grid3d.h"
#include "fe.h"
#include "fe_matr.h"
#include "poisson.h"
#include "errlog.h"

static void error_dimensions(const char *func, const char *var)
{ error_msg("%s: incorrect dimensions in %s\n", func, var); }
static void error_null(const char *func, const char *var)
{ error_msg("%s: %s is NULL \n", func, var); }
static void error_inuse(const char *func)
{ error_msg("%s: time integrator is not available.\n", func); }

static const struct fe_odens_s *glo_nsti = NULL;
static void loc_mat(const grid3d_t *grid, int el, unsigned type, double ***S);
void init_ns_time_integrator(struct fe_odens_s *nsti)
{
    if(nsti==NULL) error_null(__func__, "nsti");
    if(glo_nsti) error_inuse(__func__);
    glo_nsti = nsti;
    
    PrintLog("Begin initialization of time integrator %s.", nsti->nsode_name);
        
    InitFE(1, &nsti->elemU, nsti->grid);
    InitFE(1, &nsti->elemP, nsti->grid);
    
    int ndofsU, ndofsP;
    if(nsti->elemU == fe_bub_vec) ndofsU = fe_ndof(nsti->grid, fe_bub);
    else ndofsU = fe_ndof(nsti->grid, nsti->elemU);
    ndofsP = fe_ndof(nsti->grid, nsti->elemP);
    
    if(!nsti->u) alloc_vec(nsti->u, ndofsU);
    else if(nsti->u->N != ndofsU) error_dimensions(__func__, "u");
    if(!nsti->u0) alloc_vec(nsti->u0, ndofsU);
    else if(nsti->u0->N != ndofsU) error_dimensions(__func__, "u0");
    if(!nsti->r) alloc_vec(nsti->r, ndofsU);
    else if(nsti->r->N != ndofsU) error_dimensions(__func__, "r");
    if(!nsti->t) alloc_vec(nsti->t, ndofsU);
    else if(nsti->t->N != ndofsU) error_dimensions(__func__, "t");
    
    if(!nsti->p) alloc_scal(nsti->p, ndofsP);
    else if(nsti->p->N != ndofsP) error_dimensions(__func__, "p");
    if(!nsti->p0) alloc_scal(nsti->p0, ndofsP);
    else if(nsti->p0->N != ndofsP) error_dimensions(__func__, "p0");
    if(!nsti->phi) alloc_scal(nsti->phi, ndofsP);
    else if(nsti->phi->N != ndofsP) error_dimensions(__func__, "phi");
    if(!nsti->phi0) alloc_scal(nsti->phi0, ndofsP);
    else if(nsti->phi0->N != ndofsP) error_dimensions(__func__, "phi0");
    
    
#if defined(USE_CSC)
    nsti->mat = fe_matr_alloc(nsti->elemU, loc_mat, eo_bc_dir, nsti->grid,
            nsti->symmetric, nsti->symmetric ? "sp" : "uu");
#else
    nsti->mat = fe_matr_alloc(nsti->elemU, loc_mat, eo_bc_dir, nsti->grid,
            nsti->symmetric, NULL);
#endif

    fe_matr_structure(nsti->grid, 1, &nsti->mat);
    /* fe_matr_assemble(nsti->grid, 1, &nsti->mat); */

    nsti->poisson_solver=InitPoisson(nsti->grid, nsti->elemP);
    glo_nsti = NULL;
}

void  run_ns_time_integrator(struct fe_odens_s *nsti)
{
    
}

void done_ns_time_integrator(struct fe_odens_s *nsti)
{
    if(nsti==NULL) error_null(__func__, "nsti");
    if(glo_nsti) error_inuse(__func__);
    glo_nsti = nsti;
    
    DonePoisson(nsti->poisson_solver);
    fe_matr_free(nsti->mat);    

    PrintLog("Time integrator %s done.", nsti->nsode_name);   
}







static void loc_mat(const grid3d_t *grid, int el, unsigned type, double ***S)
{
    static double _M[MAX_LDOFS][MAX_LDOFS];
    static double *M[MAX_LDOFS];
    double **lM, **lS, **lA, **lB, **lC, **lD;    
    int i, j, nn;
    
    if(!glo_nsti) error_null(__func__, "glo_nsti");

    for(i=0;i<MAX_LDOFS;i++)
        M[i]=_M[i];
    if(S) *S=M; else return;
    
    double c[5] = {
        1,                          /* time derivative */
        glo_nsti->dt,               /* convection (a) */
        glo_nsti->dt,               /* reaction (b) */
        glo_nsti->mu*glo_nsti->dt,  /* diffusion */
        glo_nsti->dt                /* pressure */
    };
    
    if(glo_nsti->coefs)
    {
        c[0] = glo_nsti->coefs[0];                          /* time derivative */
        c[1] = glo_nsti->coefs[1]*glo_nsti->dt;             /* convection (a) */
        c[2] = glo_nsti->coefs[2]*glo_nsti->dt;             /* reaction (b) */
        c[3] = glo_nsti->coefs[3]*glo_nsti->mu*glo_nsti->dt;/* diffusion */
        c[4] = glo_nsti->coefs[4]*glo_nsti->dt;             /* pressure */
    }
    
    local_mass(grid, el, type, &lM);
    local_stif(grid, el, type, &lS);
    nn=fe_dofs(grid, el, type, NULL);
    for(i=0;i<nn;i++)
        for(j=0;j<nn;j++)
            (*S)[i][j] = c[0]*lM[i][j]+c[3]*lS[i][j];
    if(glo_nsti->a && c[1])
    {
        local_advm(grid, el, type, glo_nsti->a, &lA);
        local_advm_D(grid, el, type, glo_nsti->a, &lD);
        for(i=0;i<nn;i++)
            for(j=0;j<nn;j++)
                (*S)[i][j] += c[1]*(lA[i][j]+lD[i][j]);
    }
    if(glo_nsti->b && c[2])
    {
        local_adjoit_B(grid, el, type, glo_nsti->b, &lB);
        local_adjoit_C(grid, el, type, glo_nsti->b, &lC);
        for(i=0;i<nn;i++)
            for(j=0;j<nn;j++)
                (*S)[i][j] -= c[2]*(lB[i][j]+lC[i][j]);
    }
}




